
    GetProb = function(theta, phi){
        # theta: SNMM parameters
        # phi: nuisance parameters
        # Only applies to the case with K=1 (two time points)

        if (class(theta)=="numeric"){
          p001_p000 = theta[5]
          p010_p000 = phi[2]
          p011_p000 = theta[4] * phi[2]
          p100_p000 = theta[1] * (phi[2]*phi[5] + 1 - phi[5]) / 
            (phi[1]*phi[4] + 1 - phi[4])
          p101_p000 = theta[3] * p100_p000
          p110_p000 = phi[1] * p100_p000
          p111_p000 = theta[2] * p110_p000
          
          p_p000 = c(1, p001_p000, p010_p000, p011_p000, p100_p000, p101_p000,
                     p110_p000, p111_p000)
          p.max.index = which.max(p_p000) 
          k = p_p000 / p_p000[p.max.index] 
          g.sol = function(x){
            # prod(k) * (x^8) / prod(1-x*k) - phi[3]
            sum(log(k)) + 8 * log(x) - sum(log(1-x*k)) - log(phi[3])
          }
          # p.max = uniroot(g.sol, lower=0, upper=1, tol=1e-8)$root
          if ((g.sol(0.99))<0){
            p.max = 0.99
          }else if((g.sol(0.01))>0){
              p.max = 0.01
          }else{
            # p.max = uniroot(g.sol, lower=0, upper=1, tol=1e-8)$root #OLD
            p.max = uniroot(g.sol, lower=0.001, upper=0.999, tol=1e-8)$root #CHANGE
          }
          
          
          p = p.max * k
          names(p) = c("p000","p001","p010","p011","p100","p101","p110","p111")
          return(p)
        }else if (class(theta)=="matrix"){
          # theta <- theta.true; phi<-phi.true # DEBUG
          n.theta <- nrow(theta)
          p001_p000 = theta[,5]
          p010_p000 = phi[,2]
          p011_p000 = theta[,4] * phi[,2]
          p100_p000 = theta[,1] * (phi[,2]*phi[,5] + 1 - phi[,5]) / 
            (phi[,1]*phi[,4] + 1 - phi[,4])
          p101_p000 = theta[,3] * p100_p000
          p110_p000 = phi[,1] * p100_p000
          p111_p000 = theta[,2] * p110_p000
          
          p_p000 = cbind(rep(1,n.theta), p001_p000, p010_p000, p011_p000, p100_p000, p101_p000,
                     p110_p000, p111_p000)
          col.idx <- p.max.index <- apply(p_p000, 1,which.max)
          row.idx <- 1:n.theta
          k = p_p000 / as.numeric(p_p000[cbind(row.idx, col.idx)])
          # print(k)
          # print(phi[3])
          g.sol = function(x,ki, phi3i){
            # prod(k) * (x^8) / prod(1-x*k) - phi[3]
            # sum(log(ki)) + 8 * log(x) - sum(log(1-x*ki)) - log(phi3i[,3])
            sum(log(ki)) + 8 * log(x) - sum(log(1-x*ki)) - log(phi3i)
          }
          solve.pmax <- function(tmp){
             ki= tmp[1:8] ; phi3i = tmp[9]
             
            if ((g.sol(0.99,ki,phi3i))<0){
              return(0.99)
            }else if ((g.sol(0.01,ki,phi3i))>0){
                return(0.01)
            }else{
              # return(uniroot(g.sol,ki=ki, lower=0, upper=1, maxiter = 10)$root)
              return(uniroot(g.sol,ki=ki,phi3i=phi3i, lower=0.001, upper=0.999, tol=1e-8)$root)
              # return(uniroot(g.sol,ki=ki, lower=0.9, upper=1,extendInt ="upX", tol=1e-8)$root)
              # p.max = uniroot(g, lower=0.01, upper=1, tol=1e-8)$root
            }
          }
          
          
          
          # p.max = uniroot(g.sol, lower=0, upper=1, tol=1e-8)$root
          if (parallel.root){
            library(doParallel)
            library(doRNG)
            library(foreach)
            no_cores <- detectCores(all.tests = T) - 2
            registerDoParallel(cores=no_cores)
            p.max <-  foreach(i = 1:n.theta) %dopar% {
              ki <- k[i,]
              solve.pmax(ki)
            }
            p.max <- unlist(p.max)
          }else{
            k.and.phi3 <- cbind(k, phi[,3])
            p.max <- apply(k.and.phi3, 1, solve.pmax)
          }

          p = p.max * k
          colnames(p) = c("p000","p001","p010","p011","p100","p101","p110","p111")
          return(p)
        }
        
        
        
    }
    
    # # check 
    # theta = abs(rnorm(5,0,1))
    # phi = c(abs(rnorm(3,0,0.5)),runif(2,0,1))
    # p   = GetProb(theta,phi)
    # (phi[4] * p["p110"] + (1-phi[4]) * p["p100"]) / 
    #     (phi[5] * p["p010"] + (1-phi[5]) * p["p000"]); theta[1]
    # p["p111"] / p["p110"]; theta[2]
    # p["p101"] / p["p100"]; theta[3]
    # p["p011"] / p["p010"]; theta[4]
    # p["p001"] / p["p000"]; theta[5]
    # p["p110"] / p["p100"]; phi[1]
    # p["p010"] / p["p000"]; phi[2]
    # prod(p) / prod(1-p);   phi[3]
    
#     system.time({
#     for(i in 1:10000)    GetProb(theta,phi)
#     })
    
    

    
    
    
    